function [bv,sebv,R2v,R2vadj]=robustOLS(y,X,k,se_scheme)

%scheme: 1: Newey West (1987), 2: Hansen-Hodrick (1980)

T=size(y,1);
K=size(X,2);

bv=X\y;

e=y-X*bv;

sebv=NaN(K,1);

R2v=1-var(e)/var(y);
R2vadj=1-var(e)/var(y)*(T-1)/(T-K);

d=1/T*(X'*X);

%1: NW and 2: HH
if se_scheme == 1 || se_scheme ==2;
    
    w=e*ones(1,K).*X;
    C0=1/T*w'*w;
    S=C0;
    
    for j=1:k;
        
        Cj=1/T*w(j+1:end,:)'*w(1:end-j,:);
        
        if se_scheme ==1;
            S=S+(k+1-j)/(k+1)*(Cj+Cj');
        elseif se_scheme ==2;
            S=S+(Cj+Cj');
        end;
        
    end;
    
    varb=1/T*inv(d)*S*inv(d);
    sebv=sqrt(diag(varb));
    
end;

